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Abstract. We consider the motion of uncharged dust grains of arbitrary shape including 



cally covariant equation of motion is expressed in terms of standard optical parameters. 
Explicit expressions for secular changes of osculating orbital elements are derived in de- 
tail for the special case of the Poynting-Robertson effect. Two subcases are considered: 



5h ' (i) central acceleration due to gravity and the radial component of radiation pressure 

a ■ 
■ - ■ 

independent of the particle velocity, (ii) central acceleration given by gravity and the 
radiation force as the disturbing force. The latter case yields results which may be com- 
pared with secular orbital evolution in terms of orbital elements for an arbitrarily shaped 
dust particle. The effects of solar wind are also presented. 
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1. Introduction 

Poynting (1904) formulated a problem of finding equation of motion for a perfectly ab- 
sorbing spherical particle under the action of electromagnetic radiation. Poynting did not 
succeed in finding correct solution. The second case, closely connected with the relativis- 
tic equation of motion for a free particle under the action of electromagnetic radiation, 
was presented by Einstein (1905), who calculated the change of energy and light pres- 
sure at arbitrary angle of incidence on a plane mirror. Robertson (1937) subsequently 
derived a correct equation of motion for a perfectly absorbing spherical particle. This 
result has been applied to astronomical problems for several decades and is known as 
the Poynting- Robertson (P-R) effect. Robertson's case was relativistically generalized by 
Klacka (1992a), who showed that the generalized P-R effect holds only for the special 
case where the total momentum of the outgoing radiation per unit time is colinear with 
the incident radiation (in the proper/rest frame of reference of the particle), which may 
include radiation normally incident on a plane mirror. Since real particles interact with 
electromagnetic radiation in a complicated manner and particles of various optical prop- 
erties exist (e. g., Mishchenko et al. 2002), it is essential to have an equation of motion 
sufficiently general to cover a wide range of optical parameters, not just the limited cases 
previously investigated. As an attempt of formulating such a general equation of motion, 
we can mention Lyttleton (1976), Klacka (1994a), Klacka and Kocifaj (1994), Kocifaj et 
al. (2000). The last three presentations are of hypothetical character only, since none of 
them proves correctness of the equation of motion. Fortunately, it turned out that the 
guess is correct: Klacka (2000a) succeeded in deriving relativistically covariant equation 
of motion. This equation of motion was derived in other possible form by Klacka (2000b), 
Klacka and Kocifaj (2001), within the accuracy corresponding to Klacka (1994a). Finally, 
knowing the results obtained by Klacka and Kocifaj (1994) and Kocifaj et al. (2000) and 
having in disposal papers by Klacka (2000a, 2000b, 2000c), Kimura et al. (2002) pre- 
sented a repetition of the equation of motion corresponding to Klacka and Kocifaj (1994) 
and Kocifaj et al. (2000). 

Klacka (2000a) presented relativistically covariant form of equation of motion. Later 
on, Kimura et al. (2002) became aware that relativistic form of equation of motion is 
required to be sure that the offered equation of motion is correct. Thus, the authors 
presented relativistic derivation of the equation of motion, also. In agreement with the 
procedure presented by Klacka (2000a), the authors showed that their relativistic equation 
of motion reduces to equation of motion written in the first order in v/c (v is the velocity 
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of the particle with respect to the source of radiation and gravity). However, the process 
of completely "relativistic" considerations presented by Kimura et al. (2002) is physically 
incorrect and also the final "relativistic" result is incorrect. 

The equations of motion for a moving particle have been constructed under the as- 
sumption that thermal emission from a particle is isotropic and does not exert a radiation 
pressure force on the particle in the particle frame of reference. Recently Mishchcnko 
(2001) has formulated the radiation pressure on arbitrarily shaped particles arising from 
an anisotropy of thermal emission. In this paper, we derive the equation of motion for an 
arbitrarily shaped particle moving in a radiation field taking into account the radiation 
pressure caused by an anisotropy of thermal emission as well as scattering and absorp- 
tion of light. Relativistically covariant equation of motion is presented. Derivation of the 
equation of motion is physically fully reasoned. 

We begin by reviewing in Sec. 2 and 3 the basic physical processes in proper and 
stationary frames. The equation of motion for simultaneous action of gravity and elec- 
tromagnetic radiation is presented in Sec. 4. Sec. 5 then applies these results to the 
calculation of osculating orbital elements, including the special case of the P-R effect, 
which may be treated analytically. The calculation is carried out in Sec. 6 to first order 
in v/c and applied to the ejection of a dust particle from a parent body such as a comet 
or asteroid. Two cases are considered: (i) the disturbing acceleration is given in terms of 
velocity (Robertson 1937, Wyatt and Whipple 1950), and (ii) the electromagnetic radi- 
ation itself is a disturbing function. However, these authors (including Lyttleton 1976) 
did not obtain the correct expression for the secular change of longitude of pericenter 
(perihelion). In Sec. 7 we use the equation of motion to second order in v/c for the P-R 
effect. In particular we obtain a correct expression for the secular change of longitude of 
pericenter (perihelion). Sec. 8 then finds the secular change in the advancement of perihe- 
lion to first order in v/c, with gravity as a central acceleration. Next, we briefly discuss, 
in Sec. 9, the secular evolution of orbital elements for the P-R effect and nearly-circular 
orbits. Sec. 10 then treats the effect of the solar wind on the secular changes in the orbital 
elements to the second order in v/u, where u is the speed of solar wind particles. Finally, 
Sec. 11 summarizes our results. 

Other theoretical papers on the basic properties of the P-R effect were written during 
the last decades: Burns et al. (1979), Mediavilla and Buitrago (1989), Mignard (1992), 
Srikanth (1999), Williams (2002). Some others will be mentioned within the context of 
the discussed problems. Since some confusion exists in presenting derivations and results 
in the most cited papers (Robertson 1937; Wyatt and Whipple 1950; Burns et al. 1979; 



4 



Mignard 1992), we present detailed derivation of the secular changes of orbital elements 
for heliocentric orbits and the P-R effect. 

2. Proper reference frame of the particle stationary particle 

The term "stationary particle" will denote a particle which does not move in a given 
inertial frame of reference. Primed quantities will denote quantities measured in the 
proper reference frame of the particle - rest frame of the particle. 

The flux density of photons scattered into an elementary solid angle dQ' = 
sin 8' d6' d<j)' is proportional to p'(8',(f>') dtl 1 ', where p'(6',<p') is the "phase function". 
The phase function depends on orientation of the particle with respect to the direction 
of the incident radiation and on the particle characteristics; angles 8' , 4>' correspond to 
the direction (and orientation) of travel of the scattered radiation, 6' is the polar angle 
which vanishes for propagation along the unit vector e[ of the incident radiation. The 
phase function fulfills the normalisation condition 



The momentum of the incident beam of photons which is lost in the process of inter- 
action with the particle is proportional to the cross-section C' ext (extinction). The part 
proportional to C' abs (absorption) is emitted in the form of thermal radiation and the 
part proportional to C' ext — C' abs = C' sca is scattered. The differential scattering cross 
section dC' sca /dQ' = C' sca p'(0', <j>') depends on the polarization state of the incident light 
as well as on the incidence and scattering directions (e. g., Mishchenko et al. 2002). 

The momentum (per unit time) of the scattered photons into an elementary solid 
angle dQ' is 



S' is the flux density of radiation energy (energy flow through unit area perpendicular to 
the ray per unit time). The system of unit vectors used on the RHS of the last equation 
forms an orthogonal basis. The total momentum (per unit time) of the scattered photons 
is 




(1) 



d P ' sca = - s' c' sca P '(e',<p') k' d&v 



(2) 



where the unit vector in the direction of scattering is 



K' = cos 6' e[ + sin 6' cos 4>' e' 2 + sin 6' sin </>' e' 3 . 



(3) 




(4) 
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The momentum (per unit time) obtained by the particle due to the interaction with 
radiation - radiation force acting on the particle - is 



Y = \ S ' { C "^ < - C - £ P'(0',4>') K' dfi'} + F' e (T>) , 



(5) 



where the emission component of the radiation force acting on the particle of absolute 
temperature T" is (Mishchenko et al. 2002, pp. 63-66) 

F' e {T') = -- dJ f r' K' e (f',T >') dr' . (6) 

C JO J4n 

The unit vector r = r'/r' is given by position vector r' of the observation point with 
origin inside the particle (the emitted radiation in the far-field zone of the particle prop- 
agates in the radial direction, i. e., along the unit vector r), u)' is (angular) frequency of 
radiation, 



!>')=!£ (TV) 



K' n (r',u;')- ( Z' n (r',r",u')dr'' 



(7) 



where K' u is the (1,1) element of the particle extinction matrix, Z' u is the (1,1) element 
of the phase matrix and the Planck blackbody energy distribution is given by the well- 
known relation 

Thermal emission has to be included in the total interaction of the particle with electro- 
magnetic radiation: if the particle's absolute temperature is above zero, it can emit as 
well as scatter and absorb electromagnetic radiation. 

For the sake of brevity, we will use dimensionless efficiency factors Q' x instead of cross 
sections C' x : C' x = Q' x A' , where A' is geometrical cross section of a sphere of volume 
equal to the volume of the particle. Equation (5) can be rewritten to the form 

ITr = \ S> A> {[Q ' ext ~ < C ° S0 ' > Q ' sca] e>1 + [ ~ < S[U9 ' COS ^' > Q ' sca] e ' 2 + 

3 

[ - < sin0' sin^' > Q' sca ] e' 3 } + ^ F' ej e'j , (9) 

j=i 

where < x' >= J 47r x' p'(6' , <j>') dW and F' e - = F' e (T') ■ e!-. As for the energy, we assume 
that it is conserved: the energy (per unit time) of the incoming radiation E[, equals to 
the energy (per unit time) of the outgoing radiation (after interaction with the particle) 
E' Q . We will use the fact that time t' = r, where r is proper time. 

Summarizing important equations, we can write them in a short form 

j = l 
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where Q[ = Q' ext - < cos<9' > Q' sca , Q' 2 = - < sinfl' cos^' > Q' sca , Q' 3 = - < 
sin 9' sin</>' > Q' sca . We have added an assumption of equilibrium state when the parti- 
cle's mass does not change. 

2.1. Summary of the important equations 

Using the text concerning energy below Eq. (9) and the last Eq. (10), we may describe the 
total process of interaction in the form of the following equations (energies and momenta 
per unit time): 

E' D = E[ = A'S', 

3 

p' Q = (1 - Q[) p\ (Q' 2 e' 2 + Q' 3 e' 3 ) E'Jc ]T Kj *'j > 

Pi = (E'Jc) e[ , (11) 

The index "i" represents the incoming radiation, beam of photons, the index "o" rep- 
resents the outgoing radiation. The relation for p' Q represents a generalization of the 
following cases: 

i) Q[ = 1, Q' 2 = Q' 3 = F' el = F' e2 = F' e3 = - Robertson (1937), Robertson and Noonan 
(1968), Srikanth (1999); 

ii) Qi arbitrary, Q' 2 = Q' 3 = F' el = F' e2 = F' e3 = - Klacka (1992a); 

iii) Q[, Q' 2 , Q' 3 arbitrary, F' el = F' e2 — F' e3 = - Klacka (2000a). 

The changes of energy and momentum of the particle due to the interaction with 

electromagnetic radiation are 

d E' 
~d~r~ 
dp' 

17 



= El - E' a = , 

= P'i ~ P'o ■ (12) 



3. Stationary frame of reference 

By the term "stationary frame of reference" (laboratory frame) we mean a frame of 
reference in which particle moves with a velocity vector v = v(t). The physical quantities 
measured in the stationary frame of reference will be denoted by unprimed symbols. 

Our aim is to derive equation of motion for the particle in the stationary frame of 
reference. We will use the fact that we know this equation in the proper frame of reference 
- see Eqs. (11) and (12). 

If we have a four- vector = (A , A), where ^4° is its time component and A is its 
spatial component, generalized special Lorentz transformation yields 

A ' = 7 (A - v-A/c) , 
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A' = A + [(7 - 1) v ■ A/v 2 - 7 A°/c] v , (13) 
with inverse 

A = 7 (A ' + v ■ A'/c) , 

A = A' + [(7 - 1) v ■ A'/v 2 + 7 A°'/c] v , (14) 



where 7 = 1/ y 1 — v 2 /c 2 . 

As for four-vectors we immediately introduce four-momentum: 

p^ = (p°,p) = (E/c,p). (15) 
5.1. Incoming radiation 

Applying Eqs. (14) and (15) to quantity (E^/c^p 1 /) (four-momentum per unit time - 
proper time is a scalar quantity) and taking into account also Eqs. (11), we can write 

Ei = E[ 7 (1 + v ■ ei/c) , 

P l = ^ {e[ + [(7 - 1) v-e'Jv 2 + 7/c] v} . (16) 
Using the fact that p p = (h v, h v e\) for photons, we have 

v' = V W\ , 

e'1 = ^ {ei + [(7 - 1) v- ei /v 2 - 7/c] v} , (17) 
where 

w\ = 7 (1 — v • ei/c) . (18) 

Inserting the second of Eqs. (17) into Eqs. (18), one obtains 
E t = (1/ Wl ) E[ , 

Pl = (1/wi) (Ei/c) ei . (19) 

We have four-vector p^ = (Ei/c,Pi) = (l,ei)£',/c = (l/ioi, ei/wi) »i £i/c = 
6^ wi Ei/c. 

For monochromatic radiation the flux density of radiation energy becomes 

5" = n' hv' c ; S = nhv c , (20) 

where n and n' are concentrations of photons (photon number densities) in the corre- 
sponding reference frames. We also have continuity equation 



8 



d^r = 0, f = (cn,cne 1 ) , (21) 
with current density j M . Application of Eq. (13) then yields 

n' = wm. (22) 
Using Eqs. (17), (20) and (22) we finally obtain 

S' = w\ S . (23) 
Eqs. (11), (19) and (23) then together give E l = w 1 S A', Pi = w 1 S A' e 1 /c. 

3.2. Outgoing radiation 

The situation is analogous to that of the preceding subsection. It is only a little more 
algebraically complicated, since radiation may also spread out in directions given by unit 
vectors e2, e%. How can we find transformations for the unit vectors e' 2 and e' 3 l The 
crucial point is what physics do these unit vectors describe. The vectors e' 2 and e' 3 can 
be used to describe directions of propagation of radiation scattered by the particle. Thus, 
aberration of light also exists for each of these unit vectors. The relations between e' 2 and 
e 2 , e 3 and e 3 , are analogous to that presented by the second of Eq. (17): 

1 



e, =— { ej + [( 7 - 1) v-e./v 2 - j/c] v} , j = 1,2,3, (24) 



3 Wj 



where 

u>j=7(l - v-ej/c), j = 1,2,3. (25) 

It is worth mentioning that vectors {e'j-,j — 1,2,3} form an orthonormal set of vectors, 
and, unit vectors {ej]j = 1,2,3} arc not orthogonal unit vectors. 

Applying Eqs. (14) and (15) to the quantity (E' /c,p' ) (four-momentum per unit time 
- proper time is a scalar quantity), we can write 



E = 7 (E' a + v ■ p' ) , 



Po = Po + 



E' 

(7 - 1) v ■ p'Jv 2 + 7 -§ 



v . (26) 



Using also p\ = E[ e[/c and Eqs. (11), (24), (26), 

E = Q[ Wl Ei 7 + (1 - Qi) E % + 

wi Ei (Q' 2 + Q' 3 ) 7 - Wl Ei {Q' 2 /w 2 + Q' 3 /w 3 ) 

3 

F ej ( C M - 7 C) , 
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3 

Po = (i - Qi) — e i + Qi —^2— i v — 2^t Qj ( c e jM - 7 «) - 

.7=2 



1 3 

F;. (cej/wj-iv) . (27) 
5. 5. Equation of motion 

In analogy with Eqs. (12), we have for the changes of energy and momentum of the 
particle due to the interaction with electromagnetic radiation 

"3 — = Ei — E a , 
a t 

S? - * - * ■ m 

Putting Eqs. (27) into Eqs. (28), using also pi = (Ei/c) ei, one easily obtains 



3 = 1 



"3 



Eq. (29) may be rewritten in terms of four-vectors: 

=!£ + ;*•«)(««?-•'). <*» 

j=i \ ' 
where p^ is four-vector of the particle of mass m 

p» = mu>* , (31) 
four- vector of the world- velocity of the particle is 

/ = (7c,7«). (32) 
We have also other four-vectors 

6^ = (lM, ej M), .7 = 1,2,3. (33) 

It can be easily verified that: 

i) the quantity w Ei is a scalar quantity - see first of Eqs. (19); 

ii) Eq. (30) reduces to Eq. (10) for the case of proper inertial frame of reference of the 
particle; 

iii) Eq. (30) yields d m/d r = 0. 
We introduce 

b°j = l/wj = 7 (1 + v ■ e'j I c) , 

bj = ej I Wj = e'j + [(7 - 1) v ■ e' / v 2 + 7 / c] v , j = 1,2,3 (34) 
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for the purpose of practical calculations. Physics of these relations corresponds to aber- 
ration of light. 

We have derived an equation of motion for real dust particle under the action of 
electromagnetic radiation (including thermal emission). It is supposed that the equation 
of motion is represented by Eqs. (11) and (12) in the proper frame of reference of the 
particle. The final covariant form is represented by Eq. (30), or using Ei = w\ S A' (see 
Eqs. (11), (19) and (23)), 

j=i ^ ' 

Another form of equation of motion is presented in Appendix A. There is also ex- 
plained why considerations presented by Kimura et al. (2002) are incorrect. 
To first order in v/c, Eqs. (34)-(35) yield 

d S A! ^ 

jj = 7^7 ®'i K 1 ~ 2 v ' ei / c + v ' e J'/ c ) e i - + 

1 / r ( v • e.i\ v 
- Y F, 1 + >-) ej - - 

m ^ e M\ c I J c- 

i=i 

e, = (1 - ve'/c)e' + v/c , .7 = 1,2,3. (36) 



It is worth mentioning to stress that the values of Q'— coefficients depend on particle's 
orientation with respect to the incident radiation - their values are time dependent. A 
little different derivation of the equation of motion within the accuracy to the first order 
in v/c is presented in Klacka and Kocifaj (2001); a different set of unit vectors is used 
in Klacka and Kocifaj (1994), Kocifaj et al. (2000) and in repetition by Kimura et al. 
(2002). The set of non-orthogonal unit vectors {e^j — 1,2,3} is replaced by the set of 
orthonormal vectors {kj-,j — 1,2,3}: ei = fei, e 2 = k 2 + (v ■ k 2 /c) (k\ — k 2 ) + v/c, 
e 3 = fc 3 — (v ■ fc 3 /c) fc 3 + v/c and definition v k 3 = is used: e[ = (1 + v ■ kx/c)kx — v/c, 
e' 2 = k 2 + (v ■ k 2 /c)ki, e' 3 = k 3 . 

Papers by Klacka and Kocifaj (1994), Kocifaj et al. (2000) and Kimura et al. (2002) 
are based on the paper by Klacka (1994a), which hypothesized that quantities correspond- 
ing to components of space vector Q' = (Q'i,Q' 2 ,Q' 3 ) behave as scalars under Lorentz 
transformation. Our derivation proves correctness of the hypothesis. 

It can be verified that Eq. (35) (or Eqs. (36) within the accuracy to the first order 
in v/c) yields as special cases the situations discussed in Einstein (1905) and Robertson 
(1937) - Robertson's case is obtained simply by substituting Q[ — 1, Q' 2 — Q' 3 = 0, 
Einstein's results require a little more calculations (see Appendix B). 
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3.4- Heuristic derivation 

Since we completely understand the physics of Eq. (30), we are able to present a short 
simple derivation of Eq. (36). We have: 

i) dv'/dt = ELi {l s ' A '/( mc )} Qi + Kjl™) e 'i ( scc E q- ( 10 )); dm / dt = is supposed; 

ii) S' = 5(1 — 2v ■ ei/c) (see Eq. (23)), due to the change of concentration of photons 
n' = n(l - v ■ ei/c) - Eq. (22) and Doppler effect v' = u(l - v ■ e 1 /c) - Eq. (17); 

hi) e'j = (1 + v ■ ej/c)ej — v/c, j = 1, 2, 3 (aberration of light - see Eq. (24)). 
Taking into account these physical phenomena, we finally obtain Eq. (36). However, only 
relativistic covariant formulation proves that {Q'j , F' e - ;j = 1,2,3} behave as scalars under 
Lorentz transformation. 

A reader may compare this heuristic derivation of the general case with heuristic 
derivation of the special case (Q[ = 1, Q' 2 = Q' 3 = 0, F' e — 0) presented by Burns et al. 
(1979, pp. 5-6). 

A little different derivation of the relation between 5' and 5, based on: a) stress- 
energy tensor (energy- momentum tensor), b) transformations of electric and magnetic 
fields, may be found in Klacka (1992a: sections 2.3 and 2.5). 

3.5. Continuous distribution of density flux of energy 

For a continuous frequency distribution of density flux of energy, we can write 

Taking into account that concentration of photons fulfills n' = w\ n (Eq. (22)) and 
that Doppler effect yields v' = w\ v (Eq. (17)), we have dn' jdv' = dn/du. Lorentz 
transformation finally yields 




As a consequence, dm/dr — (this corresponds to the condition dE' /dr = 0). 
If the particle is not irradiated, then one has to use 
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instead of Eq. (38). As a consequence, Eq. (39) yields (dm/d,T) e = X^=i ^ej I c - Mass 
of the particle decreases due to the thermal emission, alone. 



4. Gravitation and electromagnetic radiation 

The generally covariant equation of motion can be immediately written on the basis of 
Eq. (38) (see e. g., Landau and Lifshitz 1975): 

where the operator D / dr is the "total" covariant derivative in the general relativistic 
sense and includes gravitational effects. 

4-1- Gravity and radiation - equation of motion to the first order in v/c 

We can write, on the basis of Eqs. (40), (25) and (33) 

d v G M G 

-J7 = T 1 ei + 

dt r z 



GM, 



3 



E*[('- ! t^) 



V 



e, - 



^ c 



+ 



\ £ K, [(l + e s - i] , (4, 



m 

3 = 1 



where 



e 3 = (1 - ve' 1 /c)e' 1 + v/c, .7 = 1,2,3, 



ft = rM / B Q {\){C' ext (\/w) - C' sca (X/w) g[(X/w)}d\ 
G M Q mc J 

^ 2 = GM Q mc I B&{X){ ~ C '^ X / W ) 92Ww)}dX , 

E>2 <-oo 

ft = ^ 1 : / B ®( A ) { - C "scaWw) 9 3 (A/ W )} dA , 



G M Q m c Jq 

w = 1 — v ■ ei / c , (42) 

if we use Sun as a source of gravitation and radiation. R Q denotes the radius of the 
Sun and Bq(X) is the solar radiance at a wavelength of A; G, M@, and r are the grav- 
itational constant, the mass of the Sun, and the distance of the particle from the cen- 
ter of the Sun, respectively. The asymmetry parameter vector g' is defined by g' = 
(l/Cg CO ) J n'(dC' sca /dCl')dCl' , where n' is a unit vector in the direction of scattering; g' 
= g[ e[ + g' 2 e' 2 + g' 3 e' 3 , e'j = (1 + v ■ ej/c) ej — v/c, e- • e' rj = 5ij. We may mention 
that f3i = 0, 02 = Q2/Q1, 03 = Q3/Q1 correspond to quantities used in Klacka and 
Kocifaj (2001). 
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5. Motion and orbital elements 

Equation of motion is given by Eq. (40), or, by Eq. (41) to the first order in v/c. We 

have to take into account that the non-dimensional parameter, for Solar System 

r^SA' Lq A' 0.02868 -, A' [m 2 ] 

G M Q mc Ql - 4nG M Q mc Ql ~ 12 it Ql m[kg] { 6) 

may change during the motion. Here L is the rate of energy outflow from the Sun, the 

solar luminosity. (3 may change and at each point of the orbit all three parameters (3\, 

f3 2 and /? 3 have to be numerically calculated (except for a very special cases, when the (3 

parameters are constant during the motion). 

If we are interested in orbital evolution in terms of osculating orbital elements, we 

come to the crucial point: Which type of osculating orbital elements is correct? Mainly, if 

the "radial radiation pressure" (the dominant term containing parameter /3 = /3i in Eq. 

(75)) has to be considered together with the central gravitational force, or not (compare 

Klacka 1992b). 

As for a definition of osculating orbital elements, we refer to any textbook on celestial 
mechanics. Brouwer and Clemence (1961) write (p. 273): "As the motion progresses under 
the influence of the various attracting bodies, the coordinates and velocity components 
at any instant may be used to obtain a set of six orbital elements. These are precisely the 
elements of the ellipse that the body would follow if from that particular instant on, the 
accelerations caused by all "perturbing" bodies ceased to exist." As for our purposes, it 
is sufficient to make a small change: "attracting bodies" are replaced by "forces" . 

In reality, our perturbing physical force corresponds to the total electromagnetic ra- 
diation force. Thus, physics incites us to use gravitational force alone, as the central ac- 
celeration. However, it may occur that the dominant term containing parameter (3\ = (3 
in Eq. (41) is comparable to central gravity term. As a consequence, orbital elements will 
change very rapidly during the motion. This would suggest to divide physical force into 
two parts and to use "radial radiation pressure" as a part of central acceleration, i. e. 
together with central gravitational acceleration. The final problem concerns the changes 
of (3 (almost random changes) which prevent us to use the term "Keplerian orbit" for an 
unperturbed problem. 

Thus, we have to use gravitational acceleration of the central body (Sun) as an ac- 
celeration defining Keplerian orbits - complete electromagnetic radiative acceleration is 
a disturbing acceleration - in Eq. (41). As for long-term orbital evolution, we may use 
mean values of orbital elements - they may be defined as time mean values of orbital 
elements when true anomaly changes in 2 it radians. 
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Lot us consider that dust particle is ejected from a parent body. Orbital elements 
of the parent body, at the moment of ejection, arc: ap, ep, ip, dp, up and 9p, i. e., 
semi-major axis, eccentricity, inclination, longitude of the ascending node, longitude of 
pcricenter - longitude of perihelion for the case of Solar System - and position angle of 
the body measured from the ascending node in the direction of the motion of the body. 
The velocity vector of the parent body, in a given coordinate system with the origin in 
the dominant central body (Sun in the case of the Solar System) is: 

vp = v R e PR + v T e PT , 

/ G Mq . 

VR = \ n 2T e P sm fp > 

y a P (1 - ep) 

/ G Mq 

vt - y aP (i - 4) (1 + ep cosfp) ' 

a P (l - ep) 
rp = - — — * — e PR 

1 + ep cos j p 

e PR = (cosfip cos9p — sin ftp sin6p cos ip , 

sin Q P cos9p + cos Op sin9p coszp , sin9p sinip) , 
epp = (— cos Op sin6p — sin Op cos9p coszp , 

— sin Op sin9p + cos Op cos9p cos ip , cos9p sinzp) , 
f P = 9 P - uop . (44) 

Radial and transversal velocity components are vr and vt, unit vector epT is normal 
to the radial vector and oriented in orientation of motion. True anomaly fp equals for 
pericenter/pcrihelion and tt for apocenter/aphelion. 

The initial conditions for the particle ejected with velocity A from the parent body 
are: 

r in = r P , 

Win = v P + A . (45) 

Eqs. (44)-(45) immediately yield that initial orbital elements of the particle are equal 
to those of the parent body for A = 0: 

a%n ap , €-%n — ep , %i n Ip , 

O in = Op , LOin = Up , Q in = 9p . (46) 

(Initial orbital elements for general case represented by Eq. (45) can be calculated from 
the relations presented in Eq. (47) below.) 
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Complete equation of motion for dust particle orbiting the central body of mass M 
under central's body electromagnetic radiation and gravity interaction is given by Eq. 
(41). We can find osculating orbital elements for dust particle according to the following 
equations: 

E= l_ v2 _GMv 

2 r 
H = r x v , 
\H\ 2 



V 



GM Q ' 

2pE 



1 + G M Q 

~2 ' 9 = a (1 - e) 



1 - e 2 
i = &rccos(H Z /\H\) , 

sinfi sini = H X /\H\ , — cosfi sini = H y /\H\ , 

e fl = r/\r\ = (x,y,z)/r , 

eAr = H/\H\ , 

e-T — e N x e R , 
sin 6 sini = z/r , 

cosQ sini = {e T ) z = (y H x - x H y ) / (r\H\) , 
sin (6 - w) = v ■ e R / (ey/G M Q /pj , 

cos (6 - u) = v ■ e T / (ey/G M Q /pj - 1 , (47) 

where q is perihelion distance. If i = 0, then 6 has to be obtained from cos(£! + 0) = x/r, 
sin(£l + 6) = y/r, and, we may take £1 in an arbitrary way (e. g., the requirement that SI 
is a continuous function of time may define for i = 0). As for the value of inclination, 
we repeat the well-known fact: prograde orbits exist for i e (0,7r/2), retrograde orbits 
exist for i e (tt/2, it). The case e > 1 describes the case when the particle is ejected from 
the Solar System. 

5.1. Poynting- Roberts on effect 

Putting Q' 2 = Q' 3 = F' el = F' e2 = F' eZ = in Eq. (40), we obtain Poynting-Robertson 
effect (Robertson 1937, Klacka 1992a): 

^ = ^Qi(^-^/c). (48) 
a t c J 

Eq. (48) enables us some analytical calculations for changes of osculating orbital 
elements: it is supposed that Q[ is a constant. We will present them, in the following 
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sections. Section 6 will consider the first order in v/c of Eq. (48), section 7 will consider 
also the second order in v/c of Eq. (48). 

6. P-R effect equation of motion to the first order in v/c 

We can write, on the basis of Eqs. (48), (41) - (43) 

where n = G Mq, e R = e.\ and (3 = /3iis a non-dimensional parameter ("the ratio of 
radiation pressure force to the gravitational force"; see also Eq. (43)). Eq. (43) reduces to 
(i — 5.7 x 10~ 5 Q[ / (g[g / cm 3 ] s[cm]), for homogeneous spherical particle: g is mass density 
and s is radius of the sphere. 

6.1. Secular changes of orbital elements - radiation pressure as a part of central acceler- 
ation 

We have to use — fi (1 — (3) eji / r 2 as a central acceleration determining osculating 
orbital elements if we want to take a time average (T is time interval between passages 
through two following pericenters) in an analytical way 



a 3 / 2 2Wo (1 ~ /?) P/j 



/3 

1 ' "* ff(//3) r 2 dfp , (50) 



a 2 , ,/l-e 2 2 ^ ->o 



^ v 

assuming non-pseudo-circular orbits and the fact that orbital elements exhibit only small 
changes during the time interval T; ap is semi-major axis, ep is eccentricity, fp is true 
anomaly, pp — ap{\ — e|); the second and the third Kepler's laws were used: r 2 dfp/dt = 
yj yu(l — (3)pp - conservation of angular momentum, ap/T 2 = fj,(l — /3)/(47r 2 ). 
Rewriting Eq. (49) into the form 



e R 



dv _ n (1 - 0) 

dt r 2 r 2 I V c J " c 

we can immediately write for components of perturbation acceleration to Keplerian mo- 
tion: 

F, fl = -2/?4^, FfiT = -pJL^Z, Fp N = 0, (52) 
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where Fp r, Fp t and Fp n are radial, transversal and normal components of perturbation 
acceleration, and two-body problem yields 



vp r = \ epsmfp , 



'0(1- 


P) 


P0 


V(i- 





vp t — y (l + e^cos/^) . (53) 

The important fact that perturbation acceleration is proportional to v/c (<C 1) ensures 
the above metioned small changes of orbital elements during the time interval T. 

Perturbation equations of celestial mechanics yield for osculating orbital elements (ap 
- semi-major axis; ep - eccentricity; ip - inclination (of the orbital plane to the reference 
frame); Clp - longitude of the ascending node; up - longitude of pericenter; Qp is the 
position angle of the particle on the orbit, when measured from the ascending node in 
the direction of the particle's motion, Qp = uop + fp): 



= i \j ..(^- n\ ^ R 6/3 S[nff3 + F f JT( - 1 + e P cos ffi)} 

dep _ I pp~ 



dt Y fi (1 - /3) 

dip r 



\FpRsinfp + Fp T L S fp + 1 t +C0S/ ; lj 
I L 1 + ep cos fp J J 







/3)P/3 




r 




dt 


vva - 


/3)P/3 


dtop 


1 / 




~df 


ep V /"(I - 
r 




VM 1 - 


/3)P/3 


dQp 


- 




dt 






where r 


= P/j/(l + 


ep cos 



w cos Qp 
sin 9^ 



Fp N 



smip 



f 17 j- it 2 + e /3 cos //3 • , I 

<^ Fp R COsfp - Fp T~j : fsin/,, f 

I 1 + C/3COS//3 J 



sin 9,3 

f> AT — — — COS ip 

smip 



yjli(\ - J) pp simp 

Inserting Eqs. (52) - (53) into Eq. (54), one easily obtains 

dap u 2ap 1 + e 2 p + 2ep cos fp + e 2 fi sin 2 fp 

~dt ~ ~ ' ^~c T~ ~ej ' 
dep /j 1 , 2 \ 

= ~ P^-{ 2e l3 + e P sm JP + lcOsfp) , 

dttp 
dt ' 

dup _ _ pjJ_ I l_ 

dt r 2 c ep 



r smQp . 

Fp N —. — — cos ip , (54) 



— = ~ P^2~ (2 - ep cos fp) sin fp , 



dQp = y/n(l-/3)pi3 
dt r 2 

It is worth mentioning that dap/dt < for any time t. 



(55) 
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The set of differential equations Eqs. (55) has to be complemented with initial con- 
ditions. If the subscript "P" denotes orbital elements of the parent body (see Eq. (44)), 
then 

_ VP 

rp in = r P = - — ■ — e PR 

1 + ep cos/p 

Vp in = v P + A , (56) 

where A is velocity with respect to the nucleus of the parent body, and 

_ P/3 in 

T p in „ Gp R in , 

1 + ep in COS fp i n 



Vp m — Vp R in ep R in + Vp T in ep T i 



in j 



P I 1 HI r 

vp p m = \ ep in smfp 



M (i 


-P) 




in 


'p (1 


-P) 


v& 


in 



V P T in = \ (1 +ep in COS fp in ) . 

y V/3 in 

ep R in = (COS ftp in COS 9/3 in - SU1 Qp in Sm Op i n COS ip in , 

sin Qp 

in COS Op in + COS Clp j n sin 6/3 i n cosip in , sin Op in s'mip ln ) , 
ep T in = {- COS VLp in sin 9/3 in - Sin 0/3 in COS 9/3 in COS ip in , 

- sin flp in sin Op m + cos flp ,„ cos 9/3 i„ cosip m , cos Op m smip in ) , 

fp in = 9/3 in ~ Up in . (57) 

We write (epM = epp x epp) 

A = Avp epp + Avp epp + Avn ep/v ■ (58) 
Using Eqs. (44), (56) - (58), we finally obtain (compare Gajdosfk and Klacka 1999): 



Vp in — i _ p 



pp (vt + Av T ) 2 + (Av N ) 2 

— 4 ' 



2 

e p in = 



f 1 + e P cos f P v^s _ \ 2 

{ i-/? 4 ) 

( 1 + e P cos/p v T s\ 2 ( vr_J_Avr\ 2 
V I-/? v T J \ v T J ' 
vt + Avt . Avn „ 

COS l p in = COSip — cos9p sin Ip , 

VTS v T s 

V T + Av T . . n 

SVCLlp in COS i Ip in = Sin Ip COS i Ip + 

VTS 

^ VN (cos9p cosip cosOp — sin9p sinOp) , 

VTS 

v T + Av T . . . 
smt p i n sin S 2/3 ^ = smip sinOp + 

VTS 

^ VN (cos9p cosip sinOp + sin9p cosOp) , 

VTS 
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ep in COS f 13 in = EE^H. (1 _|_ ep COS fp) - 1 , 

Pp 

. _ 1 + e P cos/p Vts (vr + Avr) 

Gp in Sin in 1/3 2 ' 

-L ^5 Vrp 

_ WT + Aw T . . „ 

smifl in cosB/3i„ = sinjp cosHp + 

VTS 

Av N 

cos^p , 

VTS 

sin i p in smO^m = sinip sin9p , 

v TS = {vt + Av T ) 2 + (Av N ) 2 , (59) 

where vp and vt are given by expressions presented in Eq. (78) - pp — ap (l — ep). 
For the special case A = Eq. (59) reduces to 

ap in = a P (l- 0) (l - 2p 1+ i epC ° e l fp ) " , (60) 



1 - e| - 2/3 (1 + ep cos/p) 

e0m = V — ' ( } 

where fp = Op — u>p, up i n has to be obtained from 

/? + e P cos /p 

e/3 in COS (Hp - U)p in) = — , 

in sm(6p - up in ) = - - - , (62) 

tip in = tip , iyg in =ip , Op in = @P • (63) 

Physics of Eqs. (60) - (61) is following: meteoroid escapes from the Solar System when 
the orbital energy becomes positive and this can happen when the energy due to the 
radial component of the radiation force is included, without it being necessary for this 
force to exceed the gravitational attraction (Harwit 1963). Some figures may be found 
in Klacka (1992b). As an example we may mention that particle of /3 = (1 — ep)/2 
moves in parabolic orbit if ejected at perihelion of the parent body and, in an orbit with 
eccentricity ep i„ = |1 — 3ep|/(l + ep) if released at aphelion; ep i„ = for (3 = e P = 1/3 
and aphelion ejection. 

By inserting Eqs. (55) into Eq. (50), taking into account that ep i„ < 1 (see Eq. (61)), 
one can easily obtain for the secular changes of orbital elements: 



ap 



('-4) 



dep 5 u ep 

a% [I- el) 
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As an remark we may mention that a little more simple set of differential equations than 
that represented by Eqs. (64) - (65) is obtained when semi-major axis ap is replaced by 
the quantity pp = ap ( 1 — ejg ). It can be easily verified that a new set of differential 
equations for secular evolution is given by the following equations: 

i) dpp/dt = -2/3( M /c)(l - e 2 pf' 2 /pp, 

ii) dep/dt = -(5/2)/J(/Vc)[(l - e 2 pf/ 2 /pp]ep/pp. 

These two equations immediately yield pp = pp in (ep/ep in ) 4 ^. 

The last relation corresponds to the relation between orbital angular momentum and 
eccentricity of the orbit: Hp = Hp in (epj epi n ) 2 ^ ■ 

Analogously, for secular evolution of osculating orbital energy (per unit mass) holds: 
Ep = v 2 /2 - Ml - f3)/r ee - M (l - /J)(l - e 2 )/{2pp), 
Ep = Ep in [(l - e 2 )/(l - e 2 in )] (ep in /ep)^. 
Equations i) and ii) enable to write: 

i) dpp/dt = -2/JOVc) [1 - e 2 p in (jpp/pp in fl 2 ]V 2 I pp, 

ii) dep/dt = -(5/2)/3( M /c) (efjp 2 m ) (1 - e 2 f/ 2 / ef . 

These equations point out that near-circular / pseudo-circular orbits {ep « 0) and orbits 
with small values of pp are not described by the discussed secular changes of orbital 
elements and by the considered v/c approximation for the P-R effect. 

6.2. Secular changes of orbital elements - gravitation as a central acceleration 

As it was discussed in section 5, it is not wise to use — /i (1 — (3i) en / r 2 as a central 
acceleration determining osculating orbital elements for the general case represented by 
Eqs. (40) or (41), since f}\ changes almost randomly during a motion. In order to have a 
comparable results in disposal, we have to find secular changes of orbital elements when 
central central acceleration is given by — /i e# / r 2 . Thus, we will use — /i e# / r 2 as a 
central acceleration determining osculating orbital elements. 

Taking into account Eq. (49), we take the action of electromagnetic radiation as a 
perturbation to the two-body problem. We can immediately write for components of 
perturbation acceleration: 

Fk = /?4 - 2/?4^, F T = -/?4^, F^O, (68) 
r z r A c r z c 
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where Fr, Ft and Fn are radial, transversal and normal components of perturbation 
acceleration, and two-body problem yields 



v Rd 



VTd 



'- esin/ , 
P 

(1 + ecos/) 

P 



(69) 



Perturbation equations of celestial mechanics yield for osculating orbital elements (a 
- semi-major axis; e - eccentricity; i - inclination (of the orbital plane to the reference 
frame) ; £1 - longitude of the ascending node; u> - longitude of pericenter; 6 - © = u) + / 
is the position angle of the particle on the orbit, when measured from the ascending node 
in the direction of the particle's motion): 



da 
~dt 
de 
dt 
di 
~dt 

dn 

~dt 

duo 

~dl 



dB 
dt 



2 a 



1 - e 
fp 

r 



{F R e sin/ + F T (l + e cos/)} 



/Hp 
r 

fJTp 

1 

e y 

r 

rjrp 
rjpp 



F R sin f + F T 
Fn cos 9 , 



cos/ + 



e + cos / 
1 + ecos/ 



sin G 



N 



sin i 



„ , tp 2 + ecos/ . 
Fr cos f - F t - ; ; sin / 



1 + e cos / 



sin 6 
b n —. — - COSi , 



sin % 

r 



„ sin 9 

r n —. COSl , 



where r = p/(l + e cos /). 

Inserting Eqs. (68) - (69) into Eq. (70), one easily obtains 



(70) 



da 
~dl 

de 

dt 
di 

di 
dfl 

~dt 
duj 
~dt 
dQ 



2/3 



. „ h 2o l + e 2 + 2ecos/ + e 2 sin 2 / 
2 e sin/ - fj— — 



= 



= 0, 



= 0, 



1 — e 2 r 2 c 1 

-4^ sin / - /3 4- (2e + esin 2 f + 2cos/) 



/3 



/Tip 1 



cos / — (3 4 — (2 — e cos /) sin / , 



(71) 



r 2 

It is worth mentioning that da/ dt < for any time t does not hold (perturbation corre- 
sponds to complete nongravitational acceleration, and, thus, Eqs. (55) do not hold). 

The set of differential equations Eqs. (71) has to be complemented with initial condi- 
tions. If the subscript "P" denotes orbital elements of the parent body, then the initial 
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orbital elements for the particle ejected with velocity A are: 
(v T + Av T ) 2 + (Av N ) 2 

9 ! 



Vin = VP 



V 2 

(1 + e P cos/p) - l 



COS tin 

sa\i in cosf2 in 
sini irl sinf2 in 



cos ip 



(1 + e P cos/p) 

Wp 

«t + Avt 
vp + Avt 

VTS 

Av N 



Av R 



Vt 



Av N 

VTS 

sinip cosfip + 



cosGp sinzp 



(cosBp cosip cosilp — sin6p sinfip) , 



VTS 

v T + Av T . . _ 
smip sinS2p + 

VTS 

Av N 

VTS 



cos8p cosip sin Sip + sin Op cosfip) , 



&in COS 



e in sin 
sini ira cosB 



fin = — (1 + ep cos/p) - 1 , 
VP 

/in = (1 + e P COs/p) -2 



WT + Avt . „ 

= smip cos Bp + 

Avn 

COSlp , 

VTS 

sini in sin6j n = sin ip sin6p , 

v TS = (v T + Av T ) 2 + (Av N ) 2 , 



(72) 



where vp and vt are given by expressions presented in Eq. (44) - pp — ap (l — ep). 

For the special case A = Eq. (72) reduces to a simple fact: initial osculating orbital 
elements of the particle are identical with those of the parent body: 



&in j ^in ^p , %i n Ip , 

^in — Qp ! w m = UP j @in = Qp . 



(73) 



The important fact is that Eqs. (71) contain also terms not proportional to v/c (<C 1). 
These important terms protect us to use procedure analogous to that represented by Eq. 
(50). While dispersion of osculating orbital elements is very small during a time interval 
T for the case when ji (1 — (5) is used in central acceleration, the dispersion of osculating 
orbital elements may be large during the same time interval for the case when [i is used in 
central acceleration (compare Figs. 1 and 2 in Klacka 1994b). Thus, any formal averaging 
of Eqs. (71) leading to equations analogous to Eqs. (64)-(67) is not correct. 

We have explained that it is not allowed to make a simple time averaging analogous to 
that described by Eq. (50), when — fi ep / r 2 is used as a central acceleration determining 
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osculating orbital elements. However, it is of interest if we have to numerically solve 
Newtonian vectorial equation of motion (Eq. (49)) and make numerical time averaging 
(over a time interval between passages through two following pericentcrs) , or if some 
analytical simplifications can be done - something analogous to Eqs. (64)-(67). 

Fortunately, we can make analytical calculations for the purpose of obtaining secular 
changes of semi-major axis and eccentricity even when — /i e R / r 2 is used as a central 
acceleration determining osculating orbital elements. We will derive the correct equations 
in the following two subsections. 

6.2.1. Radial forces and orbital elements 

We will proceed according to Klacka (1994), in this subsection. 
Let us consider a gravitational system of two bodies 

v = ~^e R . (74) 

Let perturbation acceleration exists in the form 

a = £ e R , (75) 
< < 1. Thus, the final equation of motion is 

Eq. (76) yields as a solution the well-known Keplerian motion and the orbit is given by 

r = ^ , (77) 

1 + e c cos (0 — uj c ) 

where 

Pc = a c (1 - ef) . (78) 

The subscript "c" denotes that orbital elements are constants of motion. If we write 

v = v cR e R + v cT e T , (79) 

where ej is a unit vector transverse to the radial vector e R in the plane of the trajectory 
(positive in the direction of motion), we have 



v c r = \J~li (1 - 0) Pc 1 e c sin (6 - u c ) , (80) 
v cT = yjfl (1 - 0) pc 1 [1 + e c cos (9 - w c )] . (81) 
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In principle, we may consider also a new set of orbital elements, which are defined by 
the central gravitational acceleration. Eqs. (77) - (81) are then of the form 

r = 1 + e cos(9 - lo) ' (82) 
p = a (1 - e 2 ) , (83) 

v = v dR e R + v dT e T , (84) 

van = P^ 1 e sin (9 - w) , (85) 



VdT = Vl^P 1 [1 + e cos (9 - lj)} ; (86) 

the fact that 9 is unchanged in both sets of orbital elements is used; e R — ( cos 9, sin 9 
), ex = ( — sin 9, cos 9 ). 

Position vector and velocity vector define a state of the body at any time. Equations 
(77) and (82) yield then 

(87) 



1 + e c cos (9 — lo c ) 1 + e cos (9 — w) 

Analogously, the other two pairs of equations (Eqs. (80) and (85), and, Eqs. (81) and 
(86)) give 



\J (I — (3) p c 1 e c sin (9 — uj c ) = \J p 1 e sin (9 — u) , 
^(1-/3) p^ 1 [1 + e c cos (9 - w c )] = Vp^ 1 [1 + e cos (9 - u>)} . (89) 
9ne can easily obtain, using Eqs. (87) and (89), 

p c (1 - 0)=p, (90) 
and Eqs. (88)-(89) yield then 

(1 - (3) e c sin (9 - lo c ) = e sin (9 - lo) , (91) 
(1 - /3) [1 + e c cos (9 - lo c )] = 1 + e cos (9 - w) . (92) 
Eq. (92) may be written in the form 

(1 - 13) e c cos (9 - oj c ) - (3 = e cos (9 - to) . (93) 
Eqs. (91) and (93) yield 

e 2 = (1 - f3) 2 e\ + (3 2 - 2 (3 (1 - (3) e c cos (9 - ia c ) . (94) 
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Equation for u is given by Eqs. (91) and (93), using also Eq. (94). Finally, Eqs. (78), 
(83), (79) and (94) yield 



f 1 + e 2 c + 2e c cos(9 - tu c )\ 
a = a c |1 + /3 1 



-l 



(95) 



Eqs. (94)-(95) show that orbital osculating elements a and e change in time, during an 
orbital revolution - the larger (3, the larger change of a and e. 

The osculating orbital elements e and a obtain values between their maxima and 
minima, which can be easily found from Eqs. (94)- (95). One can easily verify that these 
relations hold: 

e min = | (1 - /?) e c - 0\ < e < (1 - 0) e c + (3 = e max , (96) 

drain 1 <C ^ <C l^^c Q"max (Q7\ 

a c 1 - e c + (3 (1 + e c ) ~ a c ~ 1 + e c + (3 (1 - e c ) a c 
6.2.2. Mean values of semi-major axis and eccentricity 

Eqs. (96) and (97) represent interval of values for eccentricity and semi-major axis, when 
— \i en I r 2 is used as a central acceleration determining osculating orbital elements. 
However, we can make time averaging during a period T, which was decribed by Eq. 
(50): 

(e) = ^ J e(t) dt, (a) = ^ J a{t) dt , 
Eqs. (98) yield 

/\ Ci 2\3/2 1 /" 2?r e (/ c ) 



2 n Jo (1 + e c cos/ c 
<°) = - e ') 3/2 ^ ^ . a(/ :L,2 C ■ (99) 



o (1 + e c cos/ c ) 



Using Eqs. (94) and (95), we finally obtain 



/x „ 2 ,3/2 1 v/(l - /3) 2 + /?2 - 2/3 (1 - /?) e c cos f c 
{e) = ^ >*J (l + ec c s/ e ) 2 

(a) = . (1 - ej)"» f f [l + /3(l + e ; i +2eCC ° S ;i /(l ' eg)rl * • (100) 
2tt J (1 + e c cos/ c ) 

The following properties can be verified: 

i) (e) > 13, (e) = /3 for e c = 0; ii) (e)/e c > 1, (e) - e c for (3 = 0; 
iii) <9(e)/<9e c > 0; iv) 0<e)/0/3 > 0; 

v) (a) > a c /(l + (3), (a) = a c /(l + /?) for e c = 0; 

vi) (a)/a c < 1, (a) = a c for /? = 0; 

vii) d(a)/de c > 0; viii) d(a)/da c > 0; ix) d{a)/d(3 < 0. 
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6.2.3. Secular changes of semi-major axis and eccentricity 

Summarizing our results, it is possible to calculate secular evolution of eccentricity and 
semi-major axis, according to the following prescription. 
At first, initial conditions for ap and ep are calculated: 

l + epCos/p\ 1 



M in = a P (1 - /3) 1 - 2/3 



, , L 1 - e% 2(3 (1 + e P cos f P ) 

(e,) in = y 1 ^ , (101) 

supposing that particle was ejected with zero ejection velocity from a parent body - 
quantities with subscript "P" belongs to the parent body's trajectory; as for more general 
case, Eqs. (56), (58) and (59) have to be used, {ap) in = pp in / (l-e| m ), {ep) in = ep in . 

As the second step, the set of the following differential equations must be solved for 
the above presented initial conditions: 

dap _ M 2 + 3e il 

! = ,102) 
«j (i - 4) 

Finally, semi-major axis and eccentricity are calculated from: 



x ,2* [l + /? (l + e 2 p + 2epcosx)/(l-e 2 p) 

a n 1 - e p) o- / — ^ 

^ Jo (1 + e^cosx) 



da: 



, /2 i 7(1 - /3) 2 e| + /3 2 - 2/3 (1 - J) ep cos * 

e =(l- e 2)/ / V dx (103) 

^™ j (1 + e/3 cosx) 



The set of equations represented by Eqs. (101)-(103) fully corresponds to detailed 
numerical calculations of vectorial equation of motion, if we are interested in secular 
evolution of eccentricity and semi-major axis (supposing (ep)i n < 1 and ep does not 
correspond to pseudo-circular orbit) for the case when central acceleration is defined by 
gravity alone. 

It is worth mentioning that instantaneous time derivatives of semi-major axis and 
eccentricity may be both positive and negative, while secular evolution yields that semi- 
major axis and eccentricity are decreasing functions of time. 

We may mention that a little more simple procedure is obtained when semi-major 
axis ap is replaced by the quantity pp — ap ( 1 — ep ). As it was already mentioned 
behind Eq. (67), 

i) dpp/dt = -2/3( M /c) [1 - el in {pp/pp m ) 5 / 2 ] 3 / 2 / pp, 
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ii) dep/dt = -(5/2)/J(/i/c) (e S Jl/pZ m ) (1 - e 2 ) 3 / 2 / ef . 

Initial condition for pp is given by the first relation of Eq. (59), e. g., pp i n = pp / (1 — /3) 
for zero ejection velocity. The first integral of Eq. (103) is replaced by a simple relation: 
p = pp (1 - (3) (Eq. (90)). 

(It is important to stress that quantities p = (p), and a = (a), e = (e), present in Eq. 
(137), do not fulfil relation p = a(l - e 2 ).) 

7. P-R effect equation of motion to the second order in v/c 

We can write, on the basis of Eqs. (48) (Balek and Klacka 2002) 
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(104) 



where (i = G Mq, e R = e.\ and (3 = /3iis a non-dimensional parameter ("the ratio of 
radiation pressure force to the gravitational force"; see also Eq. (43)). The first term in 
Eq. (104) is Newtonian gravity for two-body problem, the second term is the Poynting- 
Robertson effect in flat spacetime, the third term corresponds to Einstein's correction to 
Newtonian gravity and the last term is a sort of interference term describing the mixing 
between the effects of gravity and radiation pressure. 



7.1. Secular changes of orbital elements - radiation pressure as a part of central acceler- 
ation 

We have to use — [i (1 — (3) e R / r 2 as a central acceleration determining osculating 
orbital elements if we want to use a fact that the elements do not change rapidly during 
a motion described by Eq. (104) - during a time interval T, where T is time interval 
between passages through two following pericenters / perihelia: 
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H -44- e R - A V -^ V -\ - 7P^^e R . (105) 
\ c / c z r \ c c J r z c A r 

We can immediately write for components of perturbation acceleration to Kcplcrian mo- 
tion, on the basis of Eq. (105) - j3 is considered to be a constant during the motion: 



"P R 



*P T 



r z c r l 2 



7 



+ 



28 



<P N 







(106) 



where Fp r, Fp t and at are radial, transversal and normal components of perturbation 
acceleration, and the two-body problem yields 
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Vp T 





P) 


Pfi 


M (1- 


P) 



pp 



ep sin fp , 

(1 + ep cos fp) , 



(107) 



where ep is osculating eccentricity of the orbit, fp is true anomaly, pp = ap(l — e|), 
is semi-major axis. The important fact that perturbation acceleration is proportional to 
v/c (<§C 1) ensures the above metioned small changes of orbital elements during the time 
interval T. 

Perturbation equations of celestial mechanics yield for osculating orbital elements (ap 
- semi-major axis; ep - eccentricity; ip - inclination (of the orbital plane to the reference 
frame); Vlp - longitude of the ascending node; up - longitude of pericenter / perihelion; 
Op is the position angle of the particle on the orbit, when measured from the ascending 
node in the direction of the particle's motion, Qp = u)p + fp): 
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sin Qp 
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r sin Qp 

—===== Fp N . . COSt/3 , 

- P)pp smip 
where r — pp/(l + ep cos fp). 

Inserting Eqs. (106) - (107) into Eq. (108), one easily obtains 
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Xa2 = ( 1 + ) e /3 sin /j8 + ( 1 + ) e | sin ( 2 //3) + 4 sin //3 



^ = -/3^i(2e /3 + e^sin 2 ^ + 2cos/ /3 ) + 



" 1 Z"' 1 -* ^xX. 1+ A- 3 



r 2 c 2 V pp \2 
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dip 
dt 
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-9 t/~ ~ f 77 x ^ul + 

H e/3 y V 2 

x w i = ^1 + cos ^ + 4e/3 ~ ( 2 ~~ 6/3 cos2 ^ + e p cos ^ 
x w2 = ^1 - ) cos ^ + 86/3 ~~ ( 6 + T~p) 6/3 cos2 ^ + ^ cos ^ 

dt r 2 

We want to find secular changes of osculating orbital elements up to 1/c 2 . 
As for the terms proportional to 1/c 2 in Eq. (109), we may take a time average (T is 
time interval between passages through two following pericenters) in an analytical way 



(109) 
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= - ] = ^ / 9(fp) r 2 dfp , (110) 

assuming non-pseudo-circular orbits and the fact that orbital elements exhibit only small 
changes during the time interval T; the second and the third Kepler's laws were used: 



r 2 dfp/dt = \J yu(l — [3)pp - conservation of angular momentum, a^/T 2 = p,(l — /3)/(47r 2 ). 
The result is: 
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Taking into account terms proportional to 1/c in Eq. (109), it is inevitable to use 
perturbation theory of the second order. Simple averaging represented by Eq. (110) is 
not sufficient: the last equation of Eq. (109) yields 
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Eq. (146) holds due to the fact, that Fp jv = - in reality, equation for dQp/dt presented 
in Eq. (108) has to be used. Instead of Eq. (110), we have more precise method of 
averaging, now. We can write, using Eq. (112): 



(g) = 9(t)dt = 



2tT 



fr-^l .if. 



271 



(113) 



x 'o V dt / j j 

We are interested in secular changes of orbital elements up to the order 1/c 2 . Within 
this accuracy we can write for the terms on the right-hand side of Eq. (113) 
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We have to find osculating orbital elements present on the right-hand sides of Eq. 
(109) in terms of true anomaly fp to the required accuracy in 1/c 2 . Differentiation of the 
relation pp = ap ( 1 — ) yields, using Eq. (109) for dap/dt and dep/dt, 



dt 



c 



|(/3 xX al + 6x X a2 ) - 2ep ^- (3 x X el + X e2 
and, using Eqs. (109), (112) and (114), we obtain 

dfp c Hp \ c Hp ep J c z pp 
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(117) 
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dep_ = _ p M 2e,3 + ep sin 2 fp + 2 cos fp 
dfp c Hp 

4 — ( I P x X el + X e2 
c 2 pp \2 



1 - P - — '- ^ sin Z/3 + 

c #3 V 
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Now, if we restrict ourselves to the first order in 1/c, Eqs. (117)-(118) reduce to 
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Eqs. (119)-(120) can be easily solved: 
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(121) 



ep = ep 



where initial values p/30, e^o correspond to /^o = 0. 

On the basis of Eq. (121), we can immediately write solution of Eq. (117): 
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On the basis of Eqs. (121) - (122), we can rewrite Eq. (118) to the form 
P VM {5 - cos2 f p) ep /2 + 2 cos fp p 2 ft 1 
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(124) 



15 3 21 

Xez = - — epo fp + - epo fp cos(2 fp) + 2fp sin fp sin fp + 

j 2 \ 2 5 

7 epo I sin(2fp) + - sm(3fp) - — ep sin(4/ /3 ) , 

4 epo J l lb 

where epo has to be inserted into expressions for X e \ and X e2 (see Eq. (109)) instead of 
ep. Solution of Eq. (124) is: 

ep = epo - -/=L= ^ -±= {2sin //3 + ^ ( 5//3 - 1 sin 2 //3 
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X inte2 = - - - ) (1 - cos/p) 
1 2 



y3 + ) e f30 (1 - cos 2/0) + Te^Q (1 - cos fp) , 



X inte s = - y e /30 /J ~ ^ e /70 (l-cos2/g) + - e [M /g sin(2/g) + 2 sin /g 

21 /7 1 \ 
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+ i (1- cos 3//,) - ^ e(30 (1- cos 4/0) . (125) 

Semi-major axis ap can be obtained from the relation pp = ap ( 1 — ej| ), using Eqs. (123) 
and (125). 

Let us calculate secular change of up to the order 1/c 2 , generated by osculating orbital 
change proportional to 1/c - compare Eqs. (109) and (111). Putting Eq. (114) into Eq. 
(113), we can write to the required accuracy, 



< * >' = U - « J * j /„ 75 I 1 + 75 If J if * ' (126 > 

since time derivative of up is proportional to 1/c, and, thus, the third term on the right- 
hand side of Eq. (148) would yield a higher order than 1/c 2 in the last integral. 
The last integral in Eq. (126) consists of two terms: 
r 2 dujp \ dup 
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Inserting relation for dup/dt from Eq. (109), we can immediately write for J w i, within 
the required accuracy: 
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where results of Eqs. (121) and (122) were used, and, Hp = /i(l — Wjpp- 
T _ a M /" 27r (2 -e^ cos/g)sin/ff _ 

7^2 — - P — / 77 «//3 = 7^21 + 7^22 , 

c Jo Hpep 

p 2n 



1.2^^ r^if^dfp 

c Jo Hp 



1.22^-2^ I" pi dfp . (129) 
c Jo H p e p 



33 



As for calculation of I U 2i, we need to consider that Hp = — [3)pp and Eq. (121) 

We obtain: 
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1 - (3 C 2 

Calculation of 1^22 is analogous to calculation of I U 2i, but we have to use Eqs. (121) and 
(122) simultaneously. The result is: 

?^^(i-±\. (131) 

1-/3 c 2 P/ 3 ^e^o e^ o y 
Eqs. (129)-(131) yield 

,„ _ iH ( ' + ' _ -l.) . (1 32) 

1 -/? c 2 p/jo y 8 2 e/30 e| 0y / 

Eqs. (127), (128) and (132) yield 

_ *L 4 ii ( i + ') . (133 ) 

We have calculated the last integral in Eq. (126) - it is represented by Eq. (133). 
Since it is proportional to 1/c 2 , the first integral in Eq. (126) has to be proportional to 
l/c°. Thus, the first integral in Eq. (126) reduces to 

Hp _ dujj£ 

r 2 dt 

Jo if/30 #/30 Jo I 1 + e /30 COS /g J 
3/2 

= "A ta (134) 

Finally, Eqs. (126), (133) and (134) lead to 

I s 3/2 

( A )J " VT^T? c 2 ^ 2 I 8 + ej • (135) 

Total shift of pericenter / perihelion is given as a sum of Eqs. (Ill) and (135). Chang- 
ing the index /30 into (3 and using pp — ap (1 — e 2 ), we finally receive: 

du^ _ 3 M 3 / 2 1 + /3 2 (-13/8 + 7/eg)/3 

( * >_ ^»/ a ri_^ (I-/?) 1 / 2 ' (136) 



It can be easily verified that < duip / dt > is 
i) an increasing function of /?, 
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ii) the perihelion circulates in a positive direction, 

iii) the rate of the advancement of perihelion is not bounded for f3 — > 1. 

Let us calculate secular change of ap to the order 1/c 2 , generated by osculating orbital 
change proportional to 1/c - compare Eqs. (109) and (111). Putting Eq. (114) into Eq. 
(113), we can write to the required accuracy, 
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since time derivatives of ap and cop are proportional to 1/c, and, thus, the third term on 
the right-hand side of Eq. (114) would yield a higher order than 1/c 2 in the last integral. 
The last integral in Eq. (137) consists of two terms: 
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Inserting relations for dap/dt and duip/dt from Eq. (109), we can immediately write for 
I a i, within the required accuracy: 
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where results of Eqs. (121) and (122) were used. 
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Considering that Hp = ^//x(l — f3)pp and Eqs. (121) and (122), we have 
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Inserting results of Eq. (141) into Eq. (140), one finally obtains 



(141) 



i a2 = - a £ 



'/30 



1 + -e^ - 



Vi - /? c vmi 



7T + (IOtt + 3) e| + (9tt + 2) e^ 
1 - e ^fO 



(142) 



We have calculated the last integral in Eq. (137) - it is represented by Eq. (142). 
Since its dominant part is proportional to 1/c, the first integral in Eq. (137) has to be 
proportional to 1/c 1 . Thus, the first integral in Eq. (137) is given, according to Eq. (115), 
as 
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It can be easily verified that the last integral in Eq. (143) equals zero, within the required 
accuracy, and 
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Finally, Eqs. (137) - (140), (142) - (144) yield 
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Total secular change of semi-major axis is given as a sum of Eqs. (Ill) and (145). 
Changing the index (30 into (3 and using Eq. (144), we can finally write 
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Let us calculate secular change of to the order 1/c 2 , generated by osculating orbital 
change proportional to 1/c - compare Eqs. (109) and (111). Putting Eq. (114) into Eq. 
(113), we can write to the required accuracy, 
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since time derivatives of ep and ujp are proportional to 1/c, and, thus, the third term on 
the right-hand side of Eq. (114) would yield a higher order than 1/c 2 in the last integral. 
The last integral in Eq. (147) consists of two terms: 
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Inserting relations for dep/dt and dujp/dt from Eq. (109), we can immediately write for 
I e i, within the required accuracy: 
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where results of Eqs. (121) and (122) were used. 



hi = - P 



f 

Jo 



2ep + ep sin 2 fp + 2 cos fp 



dfp 



(149) 



(150) 



c Jo Hp 

Considering Eq. (122) and the first relation of Eq. (141) for 1/ 'Hp, one finally obtains 
from Eq. (150): 
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Putting the results represented by Eqs. (144), (148), (149) and (151) into Eq. (147), 
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Total secular change of eccentricity is given as a sum of Eqs. (Ill) and (152). Changing 
the index 00 into in Eq. (152), we can finally write 
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Let us calculate secular change of pg (see Eq. (116)) to the order 1/c 2 , generated by 
osculating orbital change proportional to 1/c - one can easily verify that (dpg/dt)u = 0. 
Taking into account Eqs. (109) and (116), we can write to the required accuracy, 
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since time derivatives of pg and ug are proportional to 1/c, and, thus, the third term on 
the right-hand side of Eq. (114) would yield a higher order than 1/c 2 in the last integral. 
The last integral in Eq. (154) consists of two terms: 
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Inserting relations for dpg/dt and dug/dt from Eqs. (116) and (109), we can immediately 
write for I p \, within the required accuracy: 
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where Eq. (121) was used. Putting results of Eqs. (143), (144), (155) and (156) into Eq. 
(154), one obtains 
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Total secular change oipg is given as a sum of Eqs. (157) and (dpg/dt)u — 0. Changing 
the index 00 into 0, we can finally write 
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dpf> 



Comparing Eqs. (153) and (158), one easily obtains 
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This can be easily integrated and the results may be written as: 
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8. Secular change of advancement of pericenter/perihelion gravitation as a 
central acceleration 

We will use — [i e# / r 2 as a central acceleration determining osculating orbital elements. 
As we have seen in the previous subsection, corrections of the order {v/c) 2 represent very 
small corrections with respect to the order v/c- only dwp/dt is important, as for secular 
changes. Thus, the most important results, for the case when — \i e# / r 2 is used as 
a central acceleration, were presented in section 6.2. However, we have not calculated 
secular change of duj/dt in the section 6.2. We will do this now, since results of the 
preceding section 7.1 have to be used. 

We will use Eqs. (91), (93) and (94) in the form 

sin (6 — ui) = {(1 — p) e c sin (9 — u c )} / e , 
cos(0 — ui) = {(1 — p) e c cos(0 — uo c ) — p} / e , 



= yj(l - P) 



2 2 



P 2 



2/3 (1 - p) e c cos (9 - lo c ) . 



Using the fact that 6 — u = 9 — u c + u c — u = f c + uo c — cj, we can write 
dcos(9 — u) 



dt 



df c div c du> . 

— H : — > sin (9 - u>) . 

dt dt dt ' 



(161) 
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We admit that the subscript "c" may be changed into [3. Using Eq. (162), one immediately 
obtains 



du: dcos(9 - u) . i 
-dt= It >< 



duj c 
~~dt 



dfc 
dt 



(163) 



Inserting the right-hand side of the second of equations of Eq. (161) into Eq. (163), one 
finally, after differentiation, obtains 



duj doj c cos/ c 
dt dt sin f c 



_i dec 
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Differentiation of the third equation of Eq. (161) yields 
-i de 1 



dt 



ft a\ 2 dec 
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de c , dfc . , 
— cos f c + — e c sin/ c 
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Inserting Eq. (165) into Eq. (164), we can write 
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Eq. (166) is the decisive equation which enables us to find secular change of w. 

As a first approximation, let us consider Keplerian orbit characterized by conditions 
da c /dt = dp c /dt = de c /dt = duo c /dt = 0, p c = a c (l — e 2 ). Using averaging of the type of 



Eq. (98), or, Eq. (110), we can write for Eq. (166) (we use r 2 df c /dt = y/fj,(l - 0)p c ) 
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If the result ^ dx/(l + k cosx) — 2n/\/l — k 2 is used, one finally obtains 



duj 
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(168) 



^ (1 - /3)/a^ /2 if/3>e c /(l + e c ) 

Let us remind that d6 c /di = ^(1 - /3)/ac /2 (see Eq. (100)). 

Let us consider P-R effect to the first order in v/c. We will use Eq. (55) and also Eq. 
(112), which wc summarize in the following Eq. (169): 
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where Hp = y/n{l — /3)pp. Changing subscript "c" into "/3" in Eq. (166) and inserting 
the third of Eq. (169) into Eq. (166), we can write 
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On the basis of Eqs. (113)-(115), (143)-(144), (169)-(170), we can write within the re- 
quired accuracy: 
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Using also Eq. (122) in the integral of the second line of Eq. (171), we finally receive 
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(172) 



+ ecosx)' 

where i}h(x) = 1 if x > 0, "&h{x) = if x < (Heaviside's step function); it is assumed 
that (3 ^ ep/(l + ep). We see that advancement of pericenter/pcrihclion exists already in 
the first order in v/c if gravity alone is taken as a central acceleration. 



9. P-R effect and near circular orbits 

When the orbits are near circular (pseudo-circular) - central acceleration contains ra- 
diation pressure term - the orbit can not reduce in semi-major axis without increasing 
in eccentricity. Both types of orbital elements, defined by central accelerations, have 
been used in detail numerical calculations in papers by Klacka and Kaufmannova (1992, 
1993). Due to the property (e) > [3 (see section 6.2.2), one must be aware that also values 
(e) > (3 may correspond to pseudo-circular orbits. The results for pseudo-circular orbits 
were analytically confirmed by Breiter and Jackson (1998). As an advantage, the analyt- 
ical approach reproduces known results without detail numerical calculations. However, 
the analytical approach paralclly produces nonphysical results which may not be distin- 
guishable from the correct results. The nonphysical analytical results are caused by use 
of the P-R effect in the first order in v/c - very special analytical solutions will diminish 
when higher orders in v / c are used. The nonphysical analytical results have been discussed 
in more detail by Klacka (2001; see also http://xxx.lanl.gov/abs/astro-ph/0004181). 



41 



10. Solar wind effect 

We have calculated secular changes of orbital elements for the P-R effect up to the second 
order in v/c. The effects of the second order in v/c seem to be small to play an important 
role in Solar System studies. In reality similar effect coming from the Sun exists and 
this effect may play more important role. Although this effect is not connected with 
electromagnetic radiation, its secular changes of orbital elements to the first order in 
v/c (more correctly v/u, where u is the speed of solar wind particles) correspond to 
the secular changes for P-R effect. This effect is caused by solar wind particles hitting 
an interplanetary dust particle. The aim of this section is to obtain secular changes of 
orbital elements of the interplanetary dust particle under the action of the solar wind, 
up to the second order in v/u. 

Let us consider equation of motion in the form (we neglect decrease of particle's mass) 
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(173) 



where, as standardly in this paper, u = G M and non-dimensional parameter (3 is 
defined in Eq. (43) ("the ratio of radiation pressure force to the gravitational force"); 
r] w 1/3. Adding the right-hand side of Eq. (49) to the right-hand side of Eq. (173), one 
obtains 
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The first term in Eq. (174) is Newtonian gravity for two-body problem, the second term is 
the Poynting- Robertson effect in flat spacetime, the third term corresponds to the action 
of the solar wind up to the second order in v/u. 

Eq. (174) may be written in the following form: 
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10.1. Secular changes of orbital elements - radiation pressure as a part of central accel- 
eration 

Neglecting solar wind pressure term n(0/Q[)(iJ,/r 2 )(u/c), we can rewrite Eq. (175) to the 
following form: 
dv _ /i (I-/?) 

It ~ V 2 eR ~ 




J3_ jj_ I v ■ e R v 

Q[ r 2 \ c u ' 2 c u cu 



^^^{^T. + n ^-"-"-^er}, (176) 



where the second term causes deceleration of the particle's motion. We use 
— /i (1 — 0) eji I r 2 as a central acceleration determining osculating orbital elements; 
is considered to be a constant during the motion. On the basis of Eq. (176), we can 
immediately write for components of disturbing acceleration to Keplerian motion 



, , , , , V \ M vp R n vj R 



F tf_v_p? *M I vpRVp? + lizl 

V Q\ J r c Q\ r ! cu 2 c u l 
Fp n = , (177) 

where Fp r, Fp t arid Fp jv are radial, transversal and normal components of perturbation 
acceleration, and the two-body problem yields 



At (I-/?) . , 
vp r = \ epsmfp , 

V P0 



vp T = ® (1 + ep cos fp) , (178) 

V P/3 

where e/j is osculating eccentricity of the orbit, fp is true anomaly, pp = ap(l — e 2 ), ap 
is semi-major axis. The important fact that perturbation acceleration is proportional to 
v/c (<C 1) ensures small changes of orbital elements during a time interval T (T is time 
interval between passages through two following pericenters / perihelia). Thus, we may 
take a time average in an analytical way 

. (i - J) i ' : 



a A f 2tt h (1 - P) P0 

= - l r ±- f V 9(fp) r 2 dfp , (179) 

assuming non-pseudo-circular orbits and the fact that orbital elements exhibit only small 
changes during the time interval T; the second and the third Kepler's laws were used: 
r 2 dfp/dt = yj /i(l — 0)pp - conservation of angular momentum, a|/T 2 = fi(l — /3)/(47r 2 ). 



43 



Relevant perturbation equations of celestial mechanics yield for osculating orbital 
elements (ojp - longitude of pericenter / perihelion; Qp is the position angle of the particle 
on the orbit, when measured from the ascending node in the direction of the particle's 
motion, Op = u>p + fp): 
2 ap 
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M (l - (3) 
where r = pp/(l + ep cos fp). 

Putting Eq. (177) and (178) into Eq. (180), making procedure of averaging of the type 
given by Eq. (179), we finally receive 
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Initial conditions are given by Eq. (59), or, Eqs. (60) - (62). 
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Similarly, for the quantity pp = ap (l — ej^j , one can easily obtain 
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We come to the conclusion that the following relation results from Eqs. (181) and 
(182): 
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as for secular changes of pp and ep, for the simple case of equation of motion of an 
interplanetary dust particle under the action of solar wind represented by Eq. (173) and 
for the P-R effect. Thus, Eq. (182) can be solved as a separate equation 
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10.2. Secular changes of orbital elements - gravitation as a central acceleration 



We can immediately write, on the basis of Eqs. (103), (172) and (181): 
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and, for secular change of longitude of pericenter/perihelion we have 
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2tt Jo (l + ecosx) 

where z?#(x) — 1 if £ > 0, "&h{x) = if x < (Heaviside's step function); it is assumed 
that /3 ^ 6,3/(1 - e/3 ). 

Initial conditions are given by: 

i) Eq. (59), or, Eqs. (60) - (62) for the set of Eqs. (186) - (187), and, 

ii) Eq. (59), or, Eqs. (60) - (62) and Eq. (72), or, Eq. (73) for the set of Eqs. (186) and 
(188): Eqs. (72) - (73) are required for initial value of u>. 

10.3. Solar wind - discussion 

The results presented in this section hold for the most simple approximation of the solar 
wind action - only radial component of the solar wind particles is considered (Eq. (173) 
is taken as an approximation to more general equation of motion presented in Klacka 
and Saniga 1993; see also Leinert and Grim 1990). Moreover, solar wind causes decrease 
of particle's mass and the secular change of particle's mass m (present in /3, see Eq. (43)) 
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is given as dm/dt = —KA' e jf / — e|), where if is a constant depending on the 

material properties of the particle and A' e ^ is the proper effective cross sectional area of 
the particle. If A' e j ^ (area toward the Sun) is changing during the particle's motion, one 
has to use dm/dt = —KA' e ^(l — V0 R/u)/r 2 and no averaging is possible - considerations 
made in section 16 hold for the case A' e ^ = A', where A' was defined above Eq. (16), and, 
moreover, A' e ^ does not change during the particle's motion (e. g., spherical particle). 
Real velocity vector of solar wind particles is nonradial and the nonradial component 
increases with decreasing distance from the Sun (e. g., Stix 2002). As a consequence, real 
solar wind effect may cause acceleration of meteoroids in small distances from the Sun, 
instead of their deceleration. 

Finally, solar wind causes also charging of meteoroids and Lorentz force has to be 
taken into account in the case of submicron grains. 

11. Summary and conclusions 

The paper derives and presents relativistically covariant equation of motion for dust par- 
ticle under the action of electromagnetic radiation - see Eq. (40). As for most frequent 
applications to systems in the universe (e. g., meteoroids in the Solar System, dust par- 
ticles in circumstellar disks), equation of motion in the form of Eqs. (41) and (42) are 
sufficient: application of Eq. (41) (under some assumption about particle's rotation) may 
be found in Kocifaj et al. (2000), Klacka and Kocifaj (2001). Some other accelerations 
may be added to the right-hand side of Eq. (41) - e. g., gravitational perturbations of 
planets, solar wind effect (see Eq. (173) or some more precise form of equation of motion) 
or some other nongravitational accelerations. 

Special attention was devoted to the Poynting-Robertson effect, since this effect is 
standardly used in Solar System studies. We have derived secular orbital evolution for 
the P-R effect up to the second order in v/c. General equation of motion for interac- 
tion between particle and incident electromagnetic radiation shows that radiation cannot 
be considered as a part of central acceleration. The central acceleration has to contain 
only gravity of the central body (star/Sun); moreover, this corresponds to the physical 
situation when radiation effect is considered as a disturbing effect. In order to compare 
secular changes of semi- major axis and eccentricity for real particle and the P-R effect, 
the paper derives and presents also secular changes of these orbital elements for the 
P-R effect: sec Eq. (103) in section 6.2.3 and Eq. (172) in section 8 - advancement of 
pericenter/perihelion exists even in the first order of v/c. 
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Solar wind effect is also considered in section 10. Solar wind produces secular changes 
of the orbital elements analogous to the P-R effect, as for the first order in v/c - see 
first two equations in Eq. (181) and Eq. (182). However, second order in solar wind effect 
produces more significant changes of orbits than it is in the case of the P-R effect, if 
radiation pressure is a part of central acceleration - compare Eqs. (136), (146), (153) and 



Application to larger bodies, e. g., asteroids, may be found in Klacka (2000c) - some 
kind of thermal emission has to be added (quantities Q' e j, j = 1, 2, 3 present in Eqs. (40) 
- (42) have to be calculated). 

Appendix A: Another formulation of the equation of motion 

(Reference to equation of number (j) of this appendix is denoted as Eq. (A j). Reference 
to equation of number (i) of the main text is denoted as Eq. (i).) 

Proper reference frame of the particle - stationary particle 

The equation of motion of the particle in its proper frame of reference is taken in the 
form 



where E' is particle's energy, p' its momentum, r is proper time, S' is the flux density of 
radiation energy (energy flow through unit area perpendicular to the ray per unit time), 
C is the radiation pressure cross section 3x3 matrix, unit vector e\ is directed along 
the path of the incident radiation (it is supposed that beam of photons propagate in 
parallel lines) and its orientation corresponds to the orientation of light propagation; F' e 
is emission component of the radiation force acting on the particle (see Eq. (6)). 

Stationary frame of reference 

Our aim is to derive equation of motion for the particle in the stationary frame of refer- 
ence. 

Covariant equation of motion - first case 

Let the components of the pressure cross section 3x3 matrix C be an orthonormal 
basis e' bl , e' b2 , e' b3 . We may then write 



(181). 
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< = E ( e bk-e'n)e'bk , n = 1,2,3, (2) 
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where e[, e! 2 and form an ortho-normal basis. Similarly 
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On the basis of Eqs. (Al), (A2), (A3) and (A4), one obtains 
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3 = 1 

which corresponds Eq. (37). Thus, the covariant form is represented by Eq. (38). 

More straightforward consideration: Let the components of the matrix C be given in 
the orthonormal basis e[, e' 2 and e' 3 - C' kl = e k T C e[, k, I = 1 to 3. Then 

C ei = E (*? C e[) e> . (6) 

3 = 1 

Substitutions Qj A' = C e[, j = 1, 2, 3 immediately yield Eq. (37) which has already 
been rewritten to the covariant form represented by Eq. (38). 

Covariant equation of motion - second case 

We want to derive an equation of motion for the particle in the frame of reference in which 
particle moves with actual velocity v. We will use the fact that we know this equation in 
the proper frame of reference - see Eq. (Al). 

Let us have a four- vector = (A , A), where ^4° is its time component and A 
is its spatial component. Since generalized special Lorentz transformations do not form 
a group (in general, composition of two generalized special Lorentz transformations is 
not a generalized special Lorentz transformation), we will consider more general Lorentz 
transformation. This Lorentz transformation can be written as 
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A'" = A"„ A" , (7) 

where summation over repeated indices is supposed (and also in all the following equa- 
tions) - 0, 1, 2, 3 for Greek letters and 1, 2, 3 for Latin letters - and its inverse 

A" = A a " A' a . (8) 

Important property of the Lorentz transformation 

A" = A A (9) 
is, that it can be composed of the following two transformations: 
A = L A , 

A" =RA' , (10) 

where L corresponds to generalized special Lorentz transformation and R represents 
rotation in 3-dimcnsional space: 

A" = L» v A v , (11) 
and (sec Eq. (13)) 
L\ = 7, 

L\ = L\ = - 1 (v/c) l , 2 = 1,2,3, 

L'j = Sij + ( 7 - 1) (v)i (v)j /v 2 , 2 = 1,2,3, j = 1,2,3, (12) 
where Sij = 1 if i = j and Sij = if i ^ j, 

A» = LJ 1 A' a , (13) 
where (see Eq. (14)) 

'-'a 'la p II 7 ' 

T) a f3 = diag(+ 1, - 1, - 1, - 1) (14) 
and 

V = 7, 

V = V = 7(«/c)< . * = 1,2,3, 

V = % + (7 - 1) , t = 1,2,3, j = 1,2,3, (15) 
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R° — 1 , 

R\ = W = 0, * = 1,2,3 , 

W j = m , » = 1,2,3, j = 1,2,3, (16) 

where can be expressed in terms of Euler angles, and, moreover, orthogonality condi- 
tions are fulfilled: 

3 

^ ' r ij r ik = Sjk i j ' = 1, 2, 3 , k = 1, 2, 3 , 

i=l 
3 

2 r « r « =**fc ' * =1 ' 2 ' 3 ' fc = 1,2,3. (17) 

On the basis of Eqs. (A9) and (A10) we can write 
K\ = R%L\. (18) 



Eqs. (A12), (A16) and (A18) yield 
A° = £°o = 7, 

A , = L\ = - 1 {v/c) l , i = 1,2,3. 

Finally, requirement A^A^ = A ^A^ yields, 
A M A e = S g 



(19) 



(20) 



where 6$ = 1 if g = v and 5% = if g ^ v\ Eq. (A7) was also used: A^ = A^„ A" , 
4' — A "4 

Incoming radiation 

Applying Eqs. (A7) and (A19) to quantity (Ei/c,Pi) (four-momentum per unit time - 
proper time is a scalar quantity) and taking into account also the fact that pi — Ei/c e\, 
we can write 

E't = Ei Wl , (21) 
where 

ttfi = 7 (1 - v ■ ei/c) . (22) 
Using the fact that = (h v,h v ei) for photons, we have 



v' = V W\ 



(23) 
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We have four-vector = (Ei/c,Pi) = (l,e\)Ei/c = (l/tui, e\/wi) w\ Ei/c = 
bi w\ Ei/c. We have found a new four- vector b^, which is given as bf 1 = A% b\ in 
the proper frame of reference of the particle: 6/ 1 = (1, e^). The transformation of space 
components between b^ and 6/ 1 corresponds to aberration of light. 

For monochromatic radiation the flux density of radiation energy becomes 

S' = n' hv' c ; S = nhv c , (24) 

where n and n' are concentrations of photons (photon number densities) in the corre- 
sponding reference frames. We also have continuity equation 

d^.r = 0, j» = (cn,cne 1 ), (25) 
with current density Application of Eqs. (A7) and (A19) then yields 
n' = w\ n . (26) 
Using Eqs. (A23), (A24) and (A26) we finally obtain 

S' = wlS . (27) 
Eqs. (A21) and (A27) then together give E t = wi S A', p { = w x S A' ei/c. 

Covariant equation of motion 

Inspiration comes from the fact that space components of four-momentum are written 
as a product with unit vector e[ . We know that this unit vector can be generalized to a 
four-vector 

£>x = C 1 / wi,ei I wi) 

wi = 7 (1 — v ■ ei/c) . (28) 

Moreover, we know that w\ S / c is a scalar quantity - invariant of the Lorentz transfor- 
mation (see Eqs. (A27)). 

The idea is to write covariant equation of motion in the form 

d -£ = ^G^W v + \±F' ej {c^-u-) , (29) 

i=i 

if the result for the emission component of the radiation force acting on the particle was 
added (see sections 2 and 3 in the main text). The only problem is to find components 
of the tensor of the second rank G M " '. 
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In order to find G"* v , we will proceed in two steps. At first, we will rewrite Eq. (A29) 
in the proper frame of reference of the particle. Comparison with Eqs. (Al) will yield 
components of G M v . The second step is transformation from G M " to G^ v . 

In the proper frame of reference, Eq. (A29) yields 

= (ei)A + IX (ej) fc , (30) 

where the term 1 / in brackets is omitted due to the simple fact that it equals 1 in the 
proper frame of reference. Comparison with Eq. (Al) yields 

G 00 = G' j =G' j = , j = 1,2,3, 

Gf k * = -C kj , j = 1,2,3, fc= 1,2,3. (31) 

In order to find G M we have to use the Lorentz transformation 
G^ v = A Q M Kf } u G' a & , (32) 
where 

A 13 — n 7 A p 

— 'la p '/ il 7 i 

r? a/3 = diag(+ 1,- 1,- 1,- 1) . (33) 

Usage of the generalized special Lorentz transformation (Eqs. (A12) or (A15)) yields 
the results of Kimura et al. (2002) - the authors take the space tensor C as a scalar quan- 
tity under Lorentz transformation. However, generalized special Lorentz transformation 
has not to be used. The reason is that any body in torquc-free, accelerated motion un- 
dergoes rotation. This effect is known as the Thomas precession (see e. g., Robertson 
and Noonan 1968, pp. 66 - 69). Thus, process of derivation and the "relativistic" result 
presented by Kimura et al. (2002) is incorrect - they have not proved correctness of the 
euqation of motion not even in the first order in v/c. 

Consistency of the covariant formulations 

We have obtained equation of motion in the form of Eq. (A29). Another form of covariant 
equation is presented in Eq. (35) (Eq. (38)). Are these equations consistent? 

Eq. (38) is covariant equation of motion in the form corresponding to = u M /c) 
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This appendix has discussed an equation of motion of the form 

d -£) n ~ G " b ^ + it/eA^-^) ■ (35) 

We want to show that Eqs. (A34) and (A35) are equivalent, i. e., that (dp 11 / d t) 7 = 
(d^/dr) H . 

We will not write the terms F' e p in what follows. 

Multiplication of Eq. (A34) by four- vector b k M (and summation over p) yields 

^f-\ b k , = ^^ J2Q' 3 (b»b k , - ()»b kll ) , k= 1,2,3. (36) 
dT / i c j=1 

In calculations of /3 M b k ^ we will use the fact that it represents scalar product of two 
four-vectors. Thus, its value is independent on the frame of reference. For the proper 
frame of reference 

P»b kfl = l , k = 1,2,3. (37) 

It can be easily verified, that 

b k „ = 1 - e' 3 ■ e' fc = 1 - Sj k , j = 1,2,3 , k = 1,2,3 , (38) 

since in the optics of scattering processes it is assumed (defined) that unit vectors e.\ , e 2 
and e' 3 are orthogonal. Thus, we obtain (inserting results of Eqs. (A37) and (A38) into 
Eq. (A36)) 

*f) 6fc „ = - HiM Qj; , fc = 1,2,3. (39) 
a t J j c 

Multiplication of Eq. (A35) by four-vector b k M (and summation over fi) yields 
1 1 b k „ = ^G» v b lv b kfl , k = 1,2,3. (40) 



d 



T 



II 

Again, the value is independent on the frame of reference. For the proper frame of refer- 
ence Eq. (A31) yields 

,2 



a 

^ {e' k ) T (C e[) , k = 1,2,3. (41) 

It was already shown (see Eq. (A6) and the text below it) that right-hand sides of Eqs. 
(A39) and (A41) are identical. Thus, also left-hand sides of Eqs. (A39) and (A40) are 
identical: 

' <iP " " b ^= & <^> k = 1,2,3. (42) 



dr J j \d 



T 
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If we take into account that four-vectors b% may be taken in various ways, Eq. (A42) 
yields 

(£).-(£)«• 

Thus, Eqs. (A34) and (A35) are equivalent, Q. E. D.. 
Appendix B: Einstein's example 

(Reference to equation of number (j) of this appendix is denoted as Eq. (B j). Reference 
to equation of number (i) of the main text is denoted as Eq. (i).) 

Let us consider a plane mirror moving (at a given moment) along x-axis (system 
S) with velocity v = (v,Q,0), v > 0; the mirror is perpendicular to the x-axis (the 
plane of the mirror is parallel to the yz-plane). A beam of incident (hitting) photons is 
characterized by unit vector S' = (cos 9', sin 6', 0) in the proper frame (primed quantities) 
of the mirror. Reflected beam is described by the unit vector e' = (— cos 6', sin#', 0) (in 
the proper frame S'). 

The problem is: Find equation of motion of the mirror in the frame of reference S. 

Solution 1: trivial manner 

Consider one photon (frequency /') in the proper frame of the mirror. Since the directions 
(and orientations) of the incident and outgoing photons are characterized by 

S' = (+ cos 6', sin 6', 0) , 

e! = (- cos <9', sin 6>',0) , (1) 
we can immediately write 
Vt = (1,+ cos0',sin0',O) , 



p£ = (1,- cos6>',sin0',O) , (2) 

for the four-momentum of the photon before interaction with the mirror and after the 
interaction. 

As a consequence, the mirror obtains four-momentum 
p'" = vt - Vt = (0, 2 cos 9', 0, 0) . (3) 

Application of the special Lorentz transformation to Eq. (B3) yields 
P^ = ^f 2 7 (cos 9') (/3, 1,0,0) , (4) 
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where, as standardly abbreviated, 

7 = 1/ y/1 - P 2 ; (3 = v I c . (5) 

On the basis of Eq. (B4), we can immediately write equation of motion of the mirror 
dn^ R' 

^- = ^2 7 (cos 0') 09, 1,0,0) , (6) 

(IT C 

where E\ is the total energy (per unit time) of the incident radiation measured in the 
proper frame of reference. 

Solution 2: application of general theory presented in Sec. 3 of the main text 

We have to choose orthonormal vectors in the systems S': we will use e[ = S' and e' 2 
and one can easily find 

e[ = S' = (+ cos 0', sin 0', 0) , 

e' 2 = (- sin^', cos 0',O) . (7) 
We have to write (Q' 3 = 0) 

P' = ^f (Qi A + 0* 4) • (8) 
On the basis of Eqs. (B3), (B7) and (B8) we have (Q' 3 = 0) 

Qi = 2 (cos0') 2 , Q' 2 = ~ 2 (sin 61') (cos 00 . (9) 
Other prescription yields (see Eq. (34)) 
6? = 7 (l + «-e / 1 /c)=7(l + /3cos0 / ) , 

61 = ei + [(7 - 1> • e^/v 2 + 7/c] t> = (7 cos 0' + 7/?, sin 0', 0) , (10) 
6^ = 7 (1 + v ■ e' 2 /c)= 7(1-/3 sin 0') , 

62 = e' 2 + [(7 - l)t> • e' 2 /v 2 +7/c] u = (-7sin0' + 7 /3, cos 0', 0) . (11) 
Inserting Eqs. (B9) - (Bll) (and = for j = 1, 2, 3) into Eqs. (29) - (30), one obtains 

dn^ R' 

= _» |[2(cos0') 2 ] (&? -/?") + [- 2 (sin 0') (cos 0')] (&£ - /3*0} 

= M27 (cos 0')(/3, 1,0,0) . (12) 

Unit vectors e' 2 = S" and e' 2 are used. They are orthonormal in the system S'. However, 
corresponding vectors are not orthogonal in the system S. 

ex = —, {ei + [(7 - 1) v ■ e'Jv 2 + 7 /c] v} , 
w 

w' = 7 (l + t;-ei/c) , (13) 
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and analogous equation holds for vector e 2 . Inserting Eqs. (B7), one obtains. 



ei = 



1 + /3cos0''7(1 + /3cos0')' 
- sin6»' + /3 cos 6' 



cos 0' + (3 sin 9' 



} 



e 2 = 



l-/3sin<9' ' 7(1 -/?sin0')' 







} 



(14) 



It can be easily verified that scalar product of these two vectors is nonzero, in general. 
Comparison with Einstein's result 

Inserting E[ = w 2 S A' mirror cos9' into Eq. (B6) (or Eqs. (B12)) and using Eq. (B14) 
(ei = (cos 9, sin 9, 0)) for the purpose of obtaining cos 9' = (cos 9 — /3)/(l — /3 cos 9), one 
easily obtains: 

i) dE/dr — 2^S A' mirror (cos9 — fi) 2 fi; using definition of radiation pressure dE/dt = 
P v A' mirror , we have P = 2 (S / c) (cos9 - (3) 2 / (1 - /? 2 ), or, 

ii) dp/dr = 2 7 3 (S A' mirror / c) (cos9 — [3) 2 -. using definition of radiation pressure P = 
(dp/dt) I A' mirror , we have P = 2 (S / c) (cos0- pf / (1 - /3 2 ). 

Result for P is consistent with the result presented in Einstein (1905). 
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